% --------------------------------------------------------------
% Erlang module 'series'
%
% Copyright (c) Neal Binnendyk 2009, 2010.
%   <nealabq@gmail.com>
%   <http://nealabq.com/>
%
% == GPL License Notice =============================================
% This code is free software: you can redistribute it and/or modify |
% it under the terms of the GNU General Public License as published |
% by the Free Software Foundation, either version 3 of the License, |
% or (at your option) any later version.                            |
%                                                                   |
% This code is distributed in the hope that it will be useful,      |
% but WITHOUT ANY WARRANTY; without even the implied warranty of    |
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the     |
% GNU General Public License for more details. See                  |
%   <http://www.gnu.org/licenses/>                                  |
% == END License Notice =============================================
%
% This file contains a few stand-alone functions for calculating
% a numberical series. I wrote it for my blog article
%   "Infinite sums, integer sequences"
%   <http://nealabq.com/blog/2009/03/19/sums_sequences/>
%
% There are several functions here to approximate and report on
% the infinite sum for various values of A and B:
%   sum over n <- 1..infinity( (n^A)/(B^n) )
%
% For some values of A and B we know closed-form solutions for
% the infinite sum.
% --------------------------------------------------------------

-module( series).
-export( [
   report_sum_na_over_bn_numerator_for_a/1,
   report_sum_na_over_bn_numerator_for_b/1,

   list_close_sum_na_over_bn_numerator_for_a/1,
   list_close_sum_na_over_bn_numerator_for_b/1,
]).

-author( "Neal Binnendyk (nealabq)" ).
-copyright( "Copyright (c) Neal Binnendyk, Ixacta Inc, 2009, 2010" ).


% --------------------------------------------------------------
% EXPORTS - Functions to calculate terms in the series
% --------------------------------------------------------------
%   The series looks like this:
%     sum( (N^A)/(B^N) ) where N <- 1 .. infinity
%
%   A can be any number.
%   B must be > 1.
%     -1 < B =< 1 clearly does not converge.
%     B =< 0 does not work because power( X, N)
%       does not work when X <= 0 (log(0) is negative infinity).
%       Although we could change power( X, N) to
%       work when (X < 0) and N is an integer.
% --------------------------------------------------------------

% --------------------------------------------------------------
% EXPORT
% report_sum_na_over_bn_numerator_for_a( A )
%
%   For value A, report all the values for B <- 2..31
%
%   Calculates each value as the sum of the first values in an
%   infinite series.

report_sum_na_over_bn_numerator_for_a( A ) ->
  io:fwrite( "sum((n^a)/(b^n))~n"),
  io:fwrite( "  for n 1..infinity~n"),
  io:fwrite( "  for a ~w~n", [A]),
  io:fwrite( "  for b 2..31~n~n"),
  print_list(
    lists:map(
      fun( B ) -> sum_na_over_bn_numerator( A, B) end,
      lists:seq( 2, 31))).

% --------------------------------------------------------------
% EXPORT
% report_sum_na_over_bn_numerator_for_b( A )
%
%   For value B, report all the values for A <- 1..30

report_sum_na_over_bn_numerator_for_b( B ) ->
  io:fwrite( "sum((n^a)/(b^n))~n"),
  io:fwrite( "  for n 1..infinity~n"),
  io:fwrite( "  for a 1..30~n"),
  io:fwrite( "  for b ~w~n~n", [B]),
  print_list(
    lists:map(
      fun( A ) -> sum_na_over_bn_numerator( A, B) end,
      lists:seq( 1, 30))).

% --------------------------------------------------------------
% EXPORT
% list_close_sum_na_over_bn_numerator_for_a( A )
%
%   Generates the series:
%     sum( (N^A)/(B^N) ) for N <- 1..200 (infinity)
%   Multiplied by:
%     (B-1)^(A+1)
%   For values of B <- 2 .. 11 (inclusive)
%
%   For A <- 2, result is:
%     [6,12,20,30,42,56,72,90,110,132]
%   See http://www.research.att.com/~njas/sequences/A002378
%   See http://mathworld.wolfram.com/PronicNumber.html
%
%   For A <- 3, result is:
%     [26,66,132,230,366,546,776,1062,1410,1826]
%
%   For A <- 4, result is:
%     [150,480,1140,2280,4074,6720,10440,15480,22110,30624]
%
%   For A <- 1, result is:
%     [2,3,4,5,6,7,8,9,10,11]
%
%   For A <- 0, result is:
%     [1,1,1,1,1,1,1,1,1,1]
%
%   Result is not integer when A is not an integer,
%   or when A < 0.

list_close_sum_na_over_bn_numerator_for_a( A ) ->
  lists:map(
    fun( B ) -> close_sum_na_over_bn_numerator( A, B) end,
    lists:seq( 2, 31)).

% --------------------------------------------------------------
% EXPORT
% list_close_sum_na_over_bn_numerator_for_b( B )
%
%   Generates the series:
%     sum( (N^A)/(B^N) ) for N <- 1..200 (infinity)
%   Multiplied by:
%     (B-1)^(A+1)
%   For values of A <- 1 .. 10 (inclusive)
%
%   For B <- 2, result is:
%     [2,6,26,150,1082,9366,94586,1091670,14174522,204495126]
%   See http://mathworld.wolfram.com/
%         StirlingNumberoftheSecondKind.html
%   See http://www.research.att.com/~njas/sequences/A000629
%
%   For B <- 3, result is:
%     [3,12,66,480,4368,47712,608016,8855040,145083648,
%      2641216512]
%   See http://www.research.att.com/~njas/sequences/A009362
%
%   For B <- 4, result is:
%     [4,20,132,1140,12324,160020,2424132,41967540,817374564,
%      17688328020]
%
%   For B <- 5, result is:
%     [5,30,230,2280,28280,421680,7336880,145879680,3263031680,
%      81097294080]

list_close_sum_na_over_bn_numerator_for_b( B ) ->
  lists:map(
    fun( A ) -> close_sum_na_over_bn_numerator( A, B) end,
    lists:seq( 1, 30)).


% --------------------------------------------------------------
% Calculations used by exported functions above
% --------------------------------------------------------------

% --------------------------------------------------------------
% close_sum_na_over_bn_numerator( A, B )
%
%   Pretties the result from sum_na_over_bn_numerator( A, B )
%   by rounding the float numbers to nearby integers.
%   The floats have to be very close to integer values or they
%   will not come thru.

close_sum_na_over_bn_numerator( A, B ) ->
  get_close_int(
    sum_na_over_bn_numerator( A, B )).

% --------------------------------------------------------------
% sum_na_over_bn_numerator( A, B )
%
%   Returns sum( (N^A)/(B^N) ) (N <- 1..200) multiplied by
%   (B-1)^(A+1). We do this because:
%
%     sum( (N^0)/(B^N) ) is 1/((B-1)^1)
%     sum( (N^1)/(B^N) ) is B/((B-1)^2)
%     sum( (N^2)/(B^N) ) is (B*(B+1))/((B-1)^3)
%
%   We are interested in the numerator series:
%     1   B   B(B+1)  ..?..
%   The next term is NOT B(B+1)(B+2) like you might expect.
%   When B is 2 the series looks like this:
%     1   2   6   26   150   1082   9366   94586   1091670

sum_na_over_bn_numerator( A, B ) ->
  sum_na_over_bn( A, B ) * power( B - 1, A + 1).

% --------------------------------------------------------------
% sum_na_over_bn( A, B )
%
%   Sums: (N^A)/(B^N)
%     for N <- 1..200 (inclusive)
%   Tries to approximate the infinite sum by summing the first
%   200 terms.

sum_na_over_bn( A, B ) ->
  sum_a_thru_b(
    fun( N ) -> na_over_bn( A, B, N) end,
    1,
    200).

% --------------------------------------------------------------
% na_over_bn( A, B, N )
%
%   Calculates: (N^A)/(B^N)

na_over_bn( A, B, N ) ->
  power( N, A) / power( B, N).


% --------------------------------------------------------------
% Helper functions
% --------------------------------------------------------------

% --------------------------------------------------------------
% print_list( List )
%
%   Simple helper function. Prints a list of floating-point
%   numbers, but prints them as integers if they are very close
%   (within about 0.001) to an integer value.

print_list( [         ] ) -> ok;
print_list( [Float | R] ) ->
  Int = get_close_int( Float),
  io:fwrite( "~w ~n", [
    case Int of
      not_close -> Float;
      _         -> Int
    end]),
  print_list( R).


% --------------------------------------------------------------
% Math helper functions
% --------------------------------------------------------------

% --------------------------------------------------------------
% get_close_int( A )
%
%   If A is a float that is very close to an integer, returns
%   the integer. This doesn't just round, A has to be very close
%   to the integer value. Otherwise returns 'not_close'.

get_close_int( A ) when is_integer( A) -> A;
get_close_int( A ) when is_float( A) ->
  Delta = 0.001,
  Close = trunc( if A < 0 -> A - Delta; true -> A + Delta end),
  Diff  = abs( A - Close),
  case Diff =< (2 * Delta) of
    true  -> Close;
    false -> not_close
  end.

% --------------------------------------------------------------
% power( X, N )
%
%   Simple function to calculate (X^N). Uses logs like a
%   sliderule, so the answer is not exact even when A and N
%   are both integers.
%
%   Works for X > 0.
%   Works for any value of N.
%
%   Improvement: You could define this for (X == 0), as long as
%   (N /= 0) (don't forget 0.0 =/= 0).
%   This is even valid for negative numbers. ((-2)^0) == 1.
%
%   Improvement: You could define this for (X < 0) as long as
%   N is an integer. Otherwise you'd need a complex number.
%
%   Improvment: If both A and N are integers maybe this should
%   always return an integer? It's probably what people expect.
%   And the following can overflow, while overflow is unlikely
%   if we stay in integers whenever possible.

power( X, _ ) when not is_number( X) -> erlang:error( bad_arg);
power( _, N ) when not is_number( N) -> erlang:error( bad_arg);

power( X, N ) when X == 0, N == 0 -> erlang:error( 'power 0^0');
power( X, N ) when X == 0, N  < 0 -> erlang:error( 'power 1/0');
power( X, _ ) when X == 0         -> 0;

power( _, N ) when N == 0         -> 1;
power( X, N ) when N == 1         -> X;
power( X, N ) when is_integer( N) ->
  if N < 0 -> 1 / power( X, - N);
  true ->
    Small_half_n = trunc( N / 2),
    Small_sqrt   = power( X, Small_half_n),
    Small_square = Small_sqrt * Small_sqrt,
    Big_half_n   = N - Small_half_n,
    if
      Big_half_n =:=  Small_half_n      -> Small_square;
      Big_half_n =:= (Small_half_n + 1) -> Small_square * X;
    end
  end;

power( X, _ ) when X < 0 -> erlang:error( 'power X<0 complex');

power( X, N ) ->
  % ok when X > 0
  % N can be anything, although we've already handled integers
  math:exp( math:log( X) * N).

% --------------------------------------------------------------
% sum_a_thru_b( Fn, A, B )
%
%   Sums Fn( N) where N goes from A to B inclusive.
%   A and B should be integers, A <= B.
%   Fn( N) takes one integer arg and returns something
%   (probably a number) that works with lists:sum(..).
%
%   sum_a_thru_b( fun( N ) -> N end, 1, 10) returns 55
%   (which is 1+2+3+4+5+6+7+8+9+10).

sum_a_thru_b( Fn, A, B ) ->
  lists:sum( lists:map( Fn, lists:seq( A, B))).


% --------------------------------------------------------------
% Test functions
% --------------------------------------------------------------
%   Functions to calculate exact results using a known
%   closed-form solution, without summing. Closed forms for
%   when A=1 and A=2.
%
%   Calculates the formulas:
%     B/((B-1)^(A+1))
%     (B*(B+1))/((B-1)^(A+1))
% --------------------------------------------------------------
- ifndef( NODEBUG ).

-export( [
   list_exact_numerator_when_a_is_1,
   list_exact_numerator_when_a_is_2,

   exact_when_a_is_1/1,
   exact_when_a_is_2/1
]).

% --------------------------------------------------------------
% EXPORTS (debug only)
% list_exact_numerator_when_a_is_1( )
% list_exact_numerator_when_a_is_2( )
%
%   Test functions with known results.

% Returns the list:
%   [2,3,4,5,6,7,8,9,10,11]
list_exact_numerator_when_a_is_1( ) ->
  lists:map(
    fun exact_numerator_when_a_is_1/1,
    lists:seq( 2, 11)).

% Returns the list:
%   [6,12,20,30,42,56,72,90,110,132]
list_exact_numerator_when_a_is_2( ) ->
  lists:map(
    fun exact_numerator_when_a_is_2/1,
    lists:seq( 2, 11)).

% --------------------------------------------------------------
% EXPORTS (debug only)
% exact_when_a_is_1( B )
% exact_when_a_is_2( B )

exact_when_a_is_1( B ) ->
  exact_numerator_when_a_is_1( B) / exact_denominator( 1, B).

exact_when_a_is_2( B ) ->
  exact_numerator_when_a_is_2( B) / exact_denominator( 2, B).

% --------------------------------------------------------------

exact_numerator_when_a_is_1( B ) -> B.
exact_numerator_when_a_is_2( B ) -> B * (B + 1).

exact_denominator( A, B ) -> power( B - 1, A + 1).

- endif.
% --------------------------------------------------------------

% --------------------------------------------------------------
% End of file series.erl
% --------------------------------------------------------------
